Libraries
#Data Manipulation
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.2 ✓ purrr 0.3.4
## ✓ tibble 3.0.3 ✓ dplyr 1.0.1
## ✓ tidyr 1.1.1 ✓ stringr 1.4.0
## ✓ readr 1.3.1 ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
#Raster Data Handling
library(raster)
## Loading required package: sp
##
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
##
## select
## The following object is masked from 'package:tidyr':
##
## extract
#Keyless Landset Data (2013-2017)
library(getlandsat)
#Rapid Interactive Visualization
library(mapview)
## GDAL version >= 3.1.0 | setting mapviewOptions(fgb = TRUE)
#Vector Data Processing
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1
library(osmdata)
## Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
#
library(USAboundaries)
library(rnaturalearth)
library(rmapshaper)
## Registered S3 method overwritten by 'geojsonlint':
## method from
## print.location dplyr
library(readxl)
library(gghighlight)
library(ggrepel)
library(knitr)
library(ggthemes)
library(leaflet)
library(leafpop)
library(units)
## udunits system database from /Library/Frameworks/R.framework/Versions/4.0/Resources/library/units/share/udunits
Question 1:
bb = read_csv("../data/uscities.csv") %>%
filter(city == "Palo") %>%
st_as_sf(coords = c("lng", "lat"), crs = 4326) %>%
st_transform(5070) %>%
st_buffer(5000) %>%
st_bbox() %>%
st_as_sfc() %>%
st_sf()
## Parsed with column specification:
## cols(
## city = col_character(),
## city_ascii = col_character(),
## state_id = col_character(),
## state_name = col_character(),
## county_fips = col_double(),
## county_name = col_character(),
## county_fips_all = col_character(),
## county_name_all = col_character(),
## lat = col_double(),
## lng = col_double(),
## population = col_double(),
## density = col_double(),
## source = col_character(),
## military = col_logical(),
## incorporated = col_logical(),
## timezone = col_character(),
## ranking = col_double(),
## zips = col_character(),
## id = col_double()
## )
Question 2:
meta = read_csv("../data/palo-flood.csv")
## Parsed with column specification:
## cols(
## entityId = col_character(),
## acquisitionDate = col_datetime(format = ""),
## cloudCover = col_double(),
## processingLevel = col_character(),
## path = col_double(),
## row = col_double(),
## min_lat = col_double(),
## min_lon = col_double(),
## max_lat = col_double(),
## max_lon = col_double(),
## download_url = col_character()
## )
files = lsat_scene_files(meta$download_url) %>%
filter(grepl(paste0("B", 1:6, ".TIF$", collapse = "|"), file)) %>%
arrange(file) %>%
pull(file)
#"paste0" gives space between columns, "paste" will not give a space.
# We only need "6" bands because they are very large
#for (i in 1:6) download each of these bands for us
st = sapply(files, lsat_image)
## File in cache
## File in cache
## File in cache
## File in cache
## File in cache
## File in cache
# "sapply" is returning a list of vectors (instead of repeating the role 6 times)
s = stack(st) %>%
setNames(c(paste0("band", 1:6)))
# stack(st)
The dimensions of my stacked image are as follows: nrow (7811), ncol (7681), ncell (59,996,291), nlayers (6). The CRS is UTM WGS84 and the cell resolution is 30 x 30 (meters).
cropper = bb %>% st_as_sf() %>% st_transform(crs(s))
r = crop(s, cropper)
# r
The dimensions of my cropped image stacks are as follows: nrow (340), ncol (346), ncell (117,640), nlayers (6). The CRS is UTM WGS84 and the cell resolution is 30 x 30 (meters), which is the same as above.
Question 3:
#Near Infrared (NIR) wavelength is commonly used to analysis vegetation health because vegetation reflects strongly in this portion of the electromagnetic spectrum.
#Shortwave Infrared (SWIR) bands are useful for discerning what is wet and dry.
ras_name = r %>% setNames(c("coastal", "blue", "green", "red", "NIR", "SWIR 1"))
#R_G_B (natural color)
plotRGB(ras_name, r = 4, g = 3, b = 2)
#NIR_R_G (fa)(color infrared)
plotRGB(ras_name, r = 5, g = 4, b = 3)
#NIR_SWIR1_R (false color water focus)
plotRGB(ras_name, r = 5, g = 6, b = 4)
# (false color urban focus)
plotRGB(ras_name, r = 7, g = 6, b = 4)
## Warning in .local(x, ...): layer was changed to 6
par(mfrow = c(2,2))
#R_G_B (natural color)
plotRGB(ras_name, r = 4, g = 3, b = 2, stretch = "lin")
#NIR_R_G (fa)(color infrared)
plotRGB(ras_name, r = 5, g = 4, b = 3, stretch = "lin")
#NIR_SWIR1_R (false color water focus)
plotRGB(ras_name, r = 5, g = 6, b = 4, stretch = "lin")
# (false color urban focus)
plotRGB(ras_name, r = 7, g = 6, b = 4, stretch = "lin")
## Warning in .local(x, ...): layer was changed to 6
par(mfrow = c(2,2))
#R_G_B (natural color)
plotRGB(ras_name, r = 4, g = 3, b = 2, stretch = "hist")
#NIR_R_G (fa)(color infrared)
plotRGB(ras_name, r = 5, g = 4, b = 3, stretch = "hist")
#NIR_SWIR1_R (false color water focus)
plotRGB(ras_name, r = 5, g = 6, b = 4, stretch = "hist")
# (false color urban focus)
plotRGB(ras_name, r = 7, g = 6, b = 4, stretch = "hist")
## Warning in .local(x, ...): layer was changed to 6
The purpose of applying a color stretch is to seperate the variance in color to better show the information that is needed for a project. This helps with seeing and understanding patterns within a given landscape.
Question 4:
#Raster Algebra
#Normalized difference vegetation index
NDVI = (r$band5 - r$band4)/(r$band5 + r$band4)
#plot(NDVI)
#Normalized Difference Water Index
NDWI = (r$band3 - r$band5)/(r$band3 + r$band5)
#Modified, normalized difference water index
MNDWI = (r$band3 - r$band6)/(r$band3 + r$band6)
#Water ratio index
WRI = (r$band3 + r$band4)/(r$band5 + r$band6)
#Simple water index
SWI = 1/(sqrt(r$band2 - r$band6))
## Warning in sqrt(getValues(x)): NaNs produced
ras_stack = stack(NDVI, NDWI, MNDWI, WRI, SWI) %>%
setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI"))
palette = colorRampPalette(c("blue", "white", "red"))
plot(ras_stack, col = palette(256))
The images appear to deviate based on placement and quanity of blue, white, and red. For example, NDVI & SWI show the water as blue, but they differ in that SWI appears void of the color red. Another example, NDWI, MNDWI, and WRI all show the water as red, but they differ in amount of white versus blue. Additionally, all five maps have a different scale level and axis labels.
# Raster Thresholding
thres_NDVI = function(x){ifelse(x <= 0,1,0)}
thres_NDWI = function(x){ifelse(x >= 0,1,0)}
thres_MNDWI = function(x){ifelse(x >= 0,1,0)}
thres_WRI = function(x){ifelse(x >= 1,1,0)}
thres_SWI = function(x){ifelse(x <= 5,1,0)}
calc_NDVI = calc(NDVI, thres_NDVI)
calc_NDWI = calc(NDWI, thres_NDWI)
calc_MNDWI = calc(MNDWI, thres_MNDWI)
calc_WRI = calc(WRI, thres_WRI)
calc_SWI = calc(SWI, thres_SWI)
flood = stack(c(calc_NDVI, calc_NDWI, calc_MNDWI, calc_WRI, calc_SWI)) %>%
setNames(c("NDVI", "NDWI", "MNDWI", "WRI", "SWI"))
palette_fld = colorRampPalette(c("white", "blue"))
plot(flood, col = palette_fld(256))
Question 5:
set.seed(09072020)
values_g = getValues(r)
dim(values_g)
## [1] 117640 6
There is one data set (matrix?) with the size of 59,996,291 and 6 layers.
values_g = na.omit(values_g)
v5 = scale(values_g)
E5 = kmeans(v5, 12, iter.max = 100)
kmeans_raster = flood$NDVI
values(kmeans_raster) = E5$cluster
plot(kmeans_raster)
E6 = kmeans(v5, 6, iter.max = 100)
kmeans_raster_6 = flood$NDVI
clus_5 = values(kmeans_raster_6) = E6$cluster
plot(kmeans_raster_6)
E3 = kmeans(v5, 3, iter.max = 100)
kmeans_raster_3 = flood$NDVI
values(kmeans_raster_3) = E3$cluster
plot(kmeans_raster_3)
table = table(getValues(calc_MNDWI), getValues(kmeans_raster_6))
which.max(table)
## [1] 5
thres_MNDWI_5 = function(x){ifelse(x == 3,1,0)}
#calc_MNDWI = calc(MNDWI, thres_MNDWI)
calc_5 = calc(kmeans_raster_6, thres_MNDWI_5) %>% setNames(c('most flooded'))
layer_5 = addLayer(flood, calc_5)
palette_calc_5 = colorRampPalette(c("white", "blue"))
plot(layer_5, col = palette_calc_5(256))
Question 6:
total_area = cellStats(flood_5, sum)
knitr::kable(total_area, caption = "Flooded Area: Cells", col.names = c("cells"))
Flooded Area: Cells
| NDVI |
6666 |
| NDWI |
7212 |
| MNDWI |
11939 |
| WRI |
8469 |
| SWI |
15201 |
| most.flooded |
596 |
total_area_2 = total_area * (900)
knitr::kable(total_area_2, caption = "Flooded Area: Area", col.names = c("Area m2"))
Flooded Area: Area
| NDVI |
5999400 |
| NDWI |
6490800 |
| MNDWI |
10745100 |
| WRI |
7622100 |
| SWI |
13680900 |
| most.flooded |
536400 |
flood_6 = calc(layer_5, function(x){sum(x)})
plot(flood_6, col = blues9)
cell_val = ifelse(values(layer_5) ==0, NA, values(layer_5))
mapview(flood_6)
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## ellps WGS 84 in CRS definition: +proj=merc +a=6378137 +b=6378137 +lat_ts=0
## +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null +wktext +no_defs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded
## datum WGS_1984 in CRS definition
### The cell values are not an even number because they are combined to geo map data scale.
(I’m not sure why “most.flooded” was created as two: “most.flooded.1” and “most.flooded.2”, this issue gives me 7 total maps (not the 6)). There were six of the seven maps which captured the flooding in Palo, Iowa (coordinates -91.789504, 42.063058). Luckly, it appears that this was a vacant lot.
